{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Confidence Intervals (C.10)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2021-02-16T07:46:28.856351Z", "start_time": "2021-02-16T07:46:26.506270Z" } }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "from matplotlib import pyplot as plt\n", "from scipy.stats import binom\n", "from scipy import optimize" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In our article, for the sake of conciseness, we do not give the confidence interval for each number in the body of the article and in the tables." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function computes the exact confidence interval for a repeated Bernouilli, with *k* successes out of *n* independent samples. The parameter *confidence* specifies the desired level of confidence, e.g. 95%. This function is adapted from https://github.com/KazKobara/ebcic." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2021-02-16T07:46:28.865925Z", "start_time": "2021-02-16T07:46:28.857945Z" } }, "outputs": [], "source": [ "def exact_confidence_interval(n, k, confidence):\n", " \"\"\"Return exact Binomial confidence interval for the parameter.\n", "\n", " Returns:\n", " lower_p (float): lower bound of Binomial confidence interval.\n", " upper_p (float): upper bound of Binomial confidence interval.\n", " \"\"\"\n", " alpha = 1 - confidence\n", " r = alpha / 2\n", " if k > n / 2:\n", " # Cumulative error becomes large for k >> n/2,\n", " # so make k < n/2.\n", " k = n - k\n", " reverse_mode = True\n", " else:\n", " reverse_mode = False\n", "\n", " def upper(p):\n", " \"\"\"\n", " Upper bound is the p making the following 0.\n", " \"\"\"\n", " return binom.cdf(k, n, p) - r\n", "\n", " def lower(p):\n", " \"\"\"\n", " Lower bound for k>0 is the p making the following 0.\n", " \"\"\"\n", " return binom.cdf(n - k, n, 1 - p) - r\n", "\n", " if k == 0:\n", " lower_p = 0.0\n", " upper_p = 1 - (alpha)**(1 / n)\n", " else:\n", " # 0 < k <= n/2\n", " u_init = k / n\n", " l_init = k / n\n", " if k == 1:\n", " l_init = 0\n", " elif k == 2:\n", " l_init = k / (2 * n)\n", " lower_p = optimize.fsolve(lower, l_init)[0]\n", " if n == 2 and k == 1:\n", " # Exception of k/n = 1/2 and n is too small.\n", " upper_p = 1 - lower_p\n", " else:\n", " upper_p = optimize.fsolve(upper, u_init)[0]\n", "\n", " if reverse_mode:\n", " lower_p, upper_p = 1 - upper_p, 1 - lower_p\n", "\n", " return lower_p, upper_p" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For 10,000 samples, we compute the 95% confidence interval for each possible value of the number of successes *k*:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "ExecuteTime": { "end_time": "2021-02-16T07:46:52.938335Z", "start_time": "2021-02-16T07:46:28.866965Z" } }, "outputs": [], "source": [ "N_SAMPLES = 10000\n", "ks = np.array(range(N_SAMPLES + 1))\n", "lowers = []\n", "uppers = []\n", "averages = []\n", "for k in ks:\n", " lower, upper = exact_confidence_interval(n=N_SAMPLES, k=k, confidence=0.95)\n", " lowers.append(lower)\n", " uppers.append(upper)\n", " averages.append(k / N_SAMPLES)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We plot the upper and lower margins of error as a function of the empirical mean:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2021-02-16T07:46:53.170881Z", "start_time": "2021-02-16T07:46:52.939296Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhUAAAEGCAYAAADSTmfWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABb1UlEQVR4nO3dd3xUZfb48c9JDyGFhNBCCSWE3ruAQQS72BVdxbbIKq6r29T9rW31u5Zd19W1rLooKvaKKzZcIyAgvfdO6D2NJCQ5vz/uDQwxZUhmMiQ579frvmZuP/MQkjPPfYqoKsYYY4wx1RUU6ACMMcYYUzdYUmGMMcYYn7CkwhhjjDE+YUmFMcYYY3zCkgpjjDHG+ERIoAM43cXFxWmHDh0CHUadlpOTQ1RUVKDDqNOsjGuGlbP/WRn738KFC/eramJVzrWkohJNmzZlwYIFgQ6jTktPTyctLS3QYdRpVsY1w8rZ/6yM/U9Etlb1XHv8YYwxxhifsKTCGGOMMT5hSYUxxhhjfMKSCmOMMcb4RECTChE5V0TWisgGEbm3jP0iIs+6+5eJSB+PfZNEZK+IrCh1TryIfCsi693XRh777nOvtVZEzvHvpzPGGGPql4AlFSISDDwPnAd0AcaKSJdSh50HpLjLeOBFj32vA+eWcel7ge9UNQX4zl3HvfY1QFf3vBfcGIwxxhjjA5V2KRWRcOByINnzeFV9pJr3HgBsUNVN7n3eBcYAqzyOGQO8oc5UqnNFJE5EmqvqLlWdISLJZVx3DJDmvp8MpAN/dLe/q6r5wGYR2eDGMKean8MYY4wxeDdOxWfAEWAhkO/DeycB2z3WM4CBXhyTBOyq4LpNVXUXgKruEpEmHteaW8a1fkZExuPUjJCYmEh6enqFH8RUT3Z2tpWxn/myjItVyS+C/EIlrwiOFUNhsVJYfOK983ry+yKFYnWuoaoooIrz6vleIUicJVhARAj2WA8Kcl/d7aFBEBYMoUFCeDCEBgthQRAaDGFBQljwiev4m/0s+5+V8enNm6SipaqW9Zihusr6H65VOMaX93M2qr4MvAyQmpqqNtCKf9lgNv5XUsaqSk5BEUeOHuNwbgFHco9x+Ogxd/0Yh4862zLzjpGTX0RuQSE5+UXkuK+5BYXkFhQF+uOcsiCBqLAQGkaEEB0RQnREKA3DS94769HhJftDadQglEZRYSREhdEoKozo8BCvkhL7WfY/K+PTmzdJxWwR6a6qy3187wyglcd6S2BnFY4pbU/JIxIRaQ7srca1jDntFRUrh3IL2JeVz/5sd8kqYF92Pvuz8tmXnc+W3Uc5Oms6h3MLKCwuPy8PCw4irkEoMZGhRIWHEBUWTIu4MKLCg4+vNwgLoWF4CA3Cg2kQFkxESDBhIUHOEuy8hgYHEV6yzd0eEhzk1kAI4r7CyesCiDi1FYXFSlGxUqRKUZFSWFx8fL2wyNlXWFxM3rFi8o4VnXgtLOJoQRF5hcXkHyt5X0ROfhHZ+YVk5R0jO7+QQ7kFbDuYS1aesy2/sLjccgkJEhpFhRHfIIxGUaEkRIUff20SE07T6AiaxkRwKK+YwqJiQoKtY52pn7xJKoYCN4rIZpzHHwKoqvao5r3nAyki0hbYgdOI8tpSx0wFJrrtLQYCR0oebVRgKjAOeNx9/cxj+9si8jTQAqfx57xqfgZj/Kq4WNmfnc+uI3nsOnKUnYedV2c9j12Hj7InK5+iMhKFsJAgEhuG0zg6nIQIoVNyE+KjwoiNDCWuQSixkWHENSh5H0pcZBgRoUE18pigMiLOo4uaVFBYTHZ+IZlHndqbQzkFHMwp4FCu83rQY33N7kwO5hRw+OgxtFTR//aHL0loGE7TmHCaREccf01qFEnLuEiSGkXSPDaSsBBLPEzd401ScZ4/bqyqhSIyEfgaCAYmqepKEZng7n8JmAacD2wAcoGbSs4XkXdwGmQ2FpEM4EFV/Q9OMvG+iNwCbAOudK+3UkTex2kIWgjcoaq1rx7X1DmZecfYdiCXrQdy2XYwl20Hc9zXXHYfyeNY0cl/tcJDgmgeG0Hz2EgGtU+geWwETaIjaNwwnMTocBo3DKNxdPhJVfZOlXF1vwfUbWEhQcSHhBEfFeb1OceKijmQXcCezDz2ZOYxa+Fy4pq1Zk9mPnuz8th9JI9lGUc4kJN/UvIhAk3dRCPJTTSS4iJp2SiS5IQoWjaKtNoOUytVmlSo6lYR6QkMczfNVNWlvri5qk7DSRw8t73k8V6BO8o5d2w52w8AI8vZ9xjwWFXjNaaqcvIL2bgv21n25rDlQA7bD+ay9WAuh3OPnXRsQlQYreIb0LtVI5J6RB5PIJrHRtAiLpJGDUJPi9oEA6HBQTSLjaBZbAQAYfvWkJaW+rPjCgqL2XXkKDsOHSXjsPt66Cg7DueyePshpi3fddJjqZAgoVV8A5ITGpDcOIq2jaNITnBeW8RFElzDtTjGeMubLqV3Ab8EPnY3vSUiL6vqc36NzJhaRlXZl5XPhn3ZbNybzcZ9OWzY6yQSu47kHT8uOEhIioukTUIDLujenNbxDWiT0IDW8VG0io8kOiI0gJ/C+ENYSBBtEqJok1D2lN1FxcrerDy2HzzKlgM5bNnvJJ6b9+fy0+aDJzWODQsOolV8JClNounYtCEdmjqvbRtHER5iQ++YwPLm8cctwEBVzQEQkSdwxnawpMLUW3nHili7O4s1uzNZvSuL1bsyWbM7iyNHT9Q6RIUF075JQwa3S6B9k4a0T4yiQ5OGtI6Psufp5iTBQeLWRkUyoG38SftKktXNHonGpn3ZrNuTxTerdh/vphscJCQnNKBj02hSmjQkpWk0qc2iadc4yh6lmBrjTVIhgGfbgyLK7p5pTJ20NyuP5RlHWL0rk9W7nQRiy/6c47/MG4QFk9osmvO7Nye1aUM6NImmfZMomsVE2GMKU20iQpOYCJrERDCwXcJJ+/ILi9i0L4d1e7JYv8dJNNbszuLrlSeSjfCQIDo1i6ZLi1i6toiha4sYOjWLITLMajWM73mTVLwG/CQin7jrlwD/8VtExgTQ4dwClu84wrKMIyzdfpjlO46c9OiidXwDOjWL5sIeLejSPJpOzWJoHd+AIHvGbQIgPCSYzs1j6Nw85qTtececZGPtnkxW7shk5c5Mvli2k3fmbQOccTvaJTY8nmR0T4qjR8tYosK9+ZNgTPm8aaj5tIik43QtFeAmVV3s78CM8bf8wiJW7Mhk8bZDLM04wrKMw2w9kHt8f9vGUQxoG0/3pFh6tIyjc/Noa+9gaoWI0GC6tIihS4sYLu3tbFNVMg4dZeXOTFbtPMLKnZnM23yQz5Y4w/UECXRsGk2vVnH0bh1Hr1aN6NCkoTUKNaek3KRCRGJUNVNE4oEt7lKyL15VD/o/PGN851BOAQu3HmLB1kMs3HqQpRlHKHAHPEqKi6R7UixX929Fz5ZxdGsRS2wDSyBM3SHi9ChpFd+Ac7s1O779QHY+SzMOs2TbYRZvP8y05bt4d74zO0JUWDA9WsbRq3UcfVo3on9yI+IaeN/l1tQ/FdVUvA1ciDPnh2dHeXHX2/kxLmOqLeNQLnM3HWTBloMs2HqIDXuzAWduiK4tYrlhUBv6JTeiT5tGNImOCHC0xgRGQsNwzurUlLM6NQWcAdc2H8hhybbDLNnuLK/M2HS8y2tq02gGtI0/vjSNsf875oRykwpVvdB9bVtz4RhTdXuz8piz8QBzNh5g9sYDbDvoPMqIiQihX3I8l/ZOol+bRvRsFUdEqDVSM6YsQUFC+8SGtE9syOV9WwJOG42l2w8zb/NB5m05yEeLMnhz7lYA2iQ0YEByPP3bxjO4XQKt4hsEMnwTYN6MU/Gdqo6sbJsxNe1I7jHmbNrPbDeJKKmJiI4IYVC7BG46I5nB7RPo2CTaGlIaUw0RocEMbJdwvPdJYVExq3Y5bTLmbT7I9NV7+GBhBuA0Zh6a0pihHRozuF0CjU5hhFJT+1XUpiICaIAzDHYjTnQjjcGZO8OYGlVcrKzYeYQf1u4jfd0+Fm87RLE6XTr7J8dzZd+WDGnfmC4tYqxxmTF+FBIcRI+WcfRoGcetw9pRXKxs2JfNnI0HmLl+P1OX7OTtn7YhAt1axHJGh8YMS2lM3zaNrJawjquopuI24Dc4CcQij+2ZwPN+jMmY4w7mFDBz/T5+WLuPGev3sT+7ABHokRTLxBEdGN4xkZ6t4gi1wX2MCZigIKFj02g6No1m3JBkCouKWZpxhFnr9/Pjhv28OnMTL/2wkfCQIAa1S+CsTk04q1MTe1RSB1XUpuKfwD9F5E4bktvUpI37svlm5R6+WbWbJdsPowrxUWEMT2nMmamJDE9JJKFheKDDNMaUIyQ4iL5tGtG3TSPuOjuFnPxC5m0+yIz1+0hfu48Hp67kwakrSWnSkLM6NWFEpyb0bdPIvhzUAd6MdHJERG4ovVFV3/BDPKYeKlZl8bZDfLNqD9+s3M3GfTkAdE+K5a6RKaSlNqF7Uqw90jCmlooKD2GEmzw8eBFs2pfN/9bs5fu1e5n042b+PWMT0REhDO+YyKjOTTmrcxNibEyYWsmbpKK/x/sInBlAFwGWVJgqKywqZs6mA3y1YjdfLDnK4a9nExIkDGqXwLghyZzduSkt4iIDHaYxxg/aJTakXWJDbh3Wjqy8Y/y4YT//W7OX/63ZxxfLdhEaLJzRoTHndWvG2Z2bWs1kLeLNiJp3eq6LSCzwpt8iMnVWcbGyYOshPl+6k2nLd3Egp4AGYcF0aRTEL9K6MyK1iQ04ZUw9Ex0RyrndmnNut+YUFyuLtx/mqxW7+HLFbv740XKCZDkD2yZwXvdmjO7SrPILmoCqykDvuUCKrwMxdZOqsjTjCJ8v3ckXy3axOzOPiNAgRnZuykU9mpOW2oS5P84krXdSoEM1xgRYUJAcb4tx//mdWbkzk69W7ObLFbt44LOVPPDZSlLigtgevoULerQg3rqrnna8Gafic06MqBkMdAbe98XNReRc4J/udV9V1cdL7Rd3//k4ycyNqrqoonNF5D0g1b1EHHBYVXuJSDKwGljr7purqhN88TnMz207kMtHizL4ZPEOth3MJSw4iDNTE7mvRyfO7tzUJi4yxlRIROiWFEu3pFh+d04qG/Zm8eXy3bw7Zz1//mwlD3++imEpjbmkd5L9TjmNePOv8DeP94XAVlXNqO6NRSQYp2vqKCADmC8iU1V1lcdh5+HUiqQAA4EXgYEVnauqV3vc4+/AEY/rbVTVXtWN3ZQtJ7+Qact38eHCDH7afBARGNqhMXee1YHRXZsRG2mPNowxVdOhSTR3joymW1AGzTr15bMlO5m6ZAd3vbuEyNBgRnVpyiW9WzAsJdF6kQSQN20qfhCRZsAAnBqLjT669wBgg6puAhCRd4ExgGdSMQZ4Q1UVmCsicSLSHEiu7Fy3luMq4CwfxWvKUFyszN18gA8XZvDl8t0cPVZEu8ZR/P6cVC7tnWSNLY0xPiUix6d7/8M5qSzYeohPl+xg2vJdTF26k8YNw7m8TxJX9mtFhyYNAx1uvePN449bgQeA/+GMqvmciDyiqpOqee8kYLvHegZObURlxyR5ee4wYI+qrvfY1lZEFuMM4PX/VHVmWYGJyHhgPEBiYiLp6enefJ56JatAmbWjkPTtx9iTq0SGwMBmIQxNiqB9nCKSwbolGazz4lrZ2dlWxn5mZVwzrJz9r6wyHt0IzhoayvL9QczIKOSVmZv494xNdIgLYnjLEAY0CyEixLqk1wRvHn/8HuitqgcARCQBmA1UN6ko619YvTzGm3PHAu94rO8CWqvqARHpC3wqIl1VNfNnF1J9GXgZIDU1VdPS0sr+BPWMqrJw6yGm/LSNL5bvoqCwmAHJ8dw3sDXndmtW5eF309PTsTL2LyvjmmHl7H8VlfHZwN04kwt+smgH7y3YzqQVOby7rogLezTn6v6t6dM6Dqci2/iDN0lFBpDlsZ7FybUEVZUBtPJYbwns9PKYsIrOFZEQ4DKgb8k2Vc0H8t33C0VkI9ARWFDdD1LX5eQX8vGiDKb8tI01u7OIDg9hbP9WXDuwDanNogMdnjHGnKRJdAS3ndme8cPbsWjbId6bv53/LtvF+wsy6JYUww2Dkrm4Vwubh8QPKppQ7B737Q7gJxH5DKc2YAwwzwf3ng+kiEhb9x7XANeWOmYqMNFtMzEQOKKqu0RkXyXnng2s8WxQKiKJwEFVLRKRdjiNPzf54HPUWTsPH2XynC2889M2MvMK6ZYUw+OXdeeini2spbUx5rQnIvRtE0/fNvE8cFFXPlm8gzfnbOEPHy3jsWmruapfS34xqA1tEqICHWqdUdFfhpKvoBs5uXHmZ764saoWishE4GucbqGTVHWliExw978ETMPpTroBp0vpTRWd63H5azj50QfAcOARESkEioAJqnrQF5+lrlmWcZhXZ27mi+W7UFXO69acm4e2pW+bRoEOzRhjqqRheAjXD2rDLwa2Zu6mg7w5dwuTftzCq7M2k9YxkRvPaMvwlMb2aKSaKppQ7GF/31xVp+EkDp7bXvJ4r8Ad3p7rse/GMrZ9BHxUjXDrtOJi5bs1e3llxibmbTlIdHgINw1JZtyQZJtJ0BhTZ4gIg9snMLh9AruP5PH2vG28M28b4ybNo1OzaMYPb8eFPVoQFmLdUquioscfz6jqb0oNfnWcql7s18hMjSgqVqYt38Xz329gze4skuIi+fOFXbiqX0uibUIfY0wd1iw2gntGdWTiiA58tmQHr8zcxD3vL+XJr9Zy89Bkxg5obb8HT1FFjz9K5vf4WwXHmFrqWFExny7ewYvpG9m0P4cOTRryj6t7clGPFoTYwDHGmHokLCSIK/u14oq+LUlft4+Xf9jE/01bw3PfbWDswNbcOqwtTaIjAh1mrVDR44+F7siVv1TVX9RgTMaPjhUV8/6C7bzw/UZ2HD5Kl+YxvHhdH87p2owgm1rcGFOPiQgjUpswIrUJyzIO88rMzbw6cxOTZ2/h2oGtmXBme5rGWHJRkQqb8Ls9JRJFJExVC2oqKON7RcXK50t38vS369h2MJfereN49JJupKUmWsMkY4wppUfLOJ4b25vfjurI899v4I05W5ny0zbG9m/FhLT2NI+10YLL4k2/wC3AjyIyFcgp2aiqT/srKOM7qsr01Xv529drWbsni87NY3jtxv6WTBhjjBeSG0fx1JU9ufOsFF5I38CUn7bxzrztXNW/JRNHpNAs1mouPHmTVOx0lyBOdDP9WcNNc/qZt/kg/zdtNUu2H6Zt4yieG9ubC7o3t8ccxhhzilonNODxy3twx4gOvJC+kffmb+eDBRncdEZbfnVme2IbWINO8C6pWKWqH3huEJEr/RSP8YHtB3P565ermbZ8N81jI3ji8u5c3qelNcA0xphqahXfgL9e1p3b09rzj2/X8e8ZG3ln3jZuT2vPuCHJ9X6UTm+SivuAD7zYZgIsK+8Yz3+/kUmzNhMcJNwzqiO/HNaOyLD6/UNujDG+1iq+AU9f3Ytbh7Xjya/X8Ncv1/D67C3cfXZHLu/bkuB6WiNc0TgV5+GMZpkkIs967IoBCv0dmPFecbHy4cIMnvx6DfuzC7isTxJ/OKeTPeszxhg/69IihtdvGsCcjQd4/Ks1/OGjZUyes4WHLu5K/+T4QIdX4yqqqdiJM9nWxcBCj+1ZOBPBmdPA2t1Z/OmT5SzYeoi+bRrxn3H96dkqLtBhGWNMvTK4fQKf3j6E/y7bxV+nrebKl+ZwUc8W3HdeJ1rE1Z+eIhWNU7EUWCoib6vqMQARaQS0UtVDNRWgKVtuQSH//G49/5m5meiIEJ66ogdX9G1pPTqMMSZARISLerbg7M5NefGHjfz7h41MX7WHX6U5M6bWh/YW3rSp+FZELnaPXQLsE5EfVPWeik8z/vL9mr38v09XsOPwUa7q15J7z+tMfFRYoMMyxhgDRIYFc8+ojlzZtyX/N201T3+7jo8WZfDYJd0ZmtI40OH5lTfdAWJVNRO4DHhNVfviTC1uatiRo8f43QdLuen1+TQIC+a98YN48oqellAYY8xpqFV8A178RV+m3DoQAX7xn5+45/0lHMypu2NJelNTESIizYGrgD/5OR5Tjh/W7ePej5axNyufiSM6cOfIDoSH1P2qNGOMqe3O6NCYr34znH/9bwMv/bCR9LX7+POFnbmkV1Kde2TtTU3FI8DXwAZVnS8i7YD1/g3LlMjOL+S+j5cxbtI8GoaH8PGvhvC7c1ItoTDGmFokIjSY352Tyhe/HkabhAbc/d5Sxr02n91H8gIdmk9VmlSo6geq2kNVb3fXN6nq5f4PzSzPOMKFz87kvfnbmXBmez6/c6j17DDGmFostVk0H00YwsMXd2X+5oOc88wMpi7dGeiwfKbcpEJE/uC+Piciz5ZefHFzETlXRNaKyAYRubeM/eLeb4OILBORPpWdKyIPicgOEVniLud77LvPPX6tiJzji8/gD8XFyqszN3HZiz+SX1jMu+MHc+95nepFy2FjjKnrgoKEcUOSmXbXMNolRvHrdxYz8e1FHM6t/W0tKmpTsdp9XeCPG7vTqj8PjAIygPkiMlVVV3kcdh6Q4i4DgReBgV6c+w9V/Vup+3UBrgG6Ai2A6SLSUVWL/PH5qupAdj6/+2Ap36/dx+guTXnyih7ENbCGmMYYU9e0bRzFB7cN5qUfNvLM9PXM23yQp6/qVat7iFQ0TsXn7utkP917AE47jU0AIvIuMAbwTCrGAG+oqgJzRSTObTSa7MW5pY0B3lXVfGCziGxwY5jj249VdcsyDnPbmws5kFPAI2O6cv2gNnWuEY8xxpgTQoKDmHhWCmmpTfjNe0u4ftJP3DmiA3ed3bFWDvVd0TDdn1PBbKSqenE1750EbPdYz8CpjajsmCQvzp0oIjfg1LL81h2sKwmYW8a1fkZExgPjARITE0lPT/fuE1XDrB3HeH1lAbFhwp8GhNM6fws//LDF7/c9HWRnZ9dIGddnVsY1w8rZ/+pyGf++h/Lm6hCe/d8GvlmyiQk9womLqF0TQVb0+KPk8cFlQDPgLXd9LLDFB/cuKwUrncSUd0xF574I/MVd/wvwd+BmL+/nbFR9GXgZIDU1VdPS0so6zCeOFRXz2BereX35Fs7okMBzY/vUu3En0tPT8WcZGyvjmmLl7H91vYzPORs+XJjBnz9dwaMLinjm6h616nFIRY8/fgAQkb+o6nCPXZ+LyAwf3DsDaOWx3hJnvhFvjgkr71xV3VOyUUReAf57CverUUeOHmPCmwuZs+kAvxzWlj+e28mmJzfGmHruir4t6dkyltunLOKGST9x//mduWVo21rxONybv2CJ7tgUAIhIWyDRB/eeD6SISFsRCcNpRDm11DFTgRvcXiCDgCOququic902FyUuBVZ4XOsaEQl3P0MKMM8Hn6NKMg7lcsWLs1mw9SBPX9WTP13QxRIKY4wxAKQ0jebTO85gVJemPPrFan7/4TLyC0+rfgVl8mZEzbuBdBHZ5K4n47Y3qA5VLRSRiTgDawUDk1R1pYhMcPe/BEzDmX59A5AL3FTRue6lnxSRXjiPNrYAt7nnrBSR93EacxYCdwSq58eKHUe46fX55B0rYvLNAxjSvvZUbRljjKkZUeEhvHhdX/753Xr++d16Nu7L5t+/6EuTmIhAh1auSpMKVf1KRFKATu6mNW4PimpT1Wk4iYPntpc83itwh7fnutuvr+B+jwGPVTVeX/hxw35++cYCGjUIY8qtA+nYNDqQ4RhjjDmNBQUJd4/qSGqzaH77/lIuef5HJt88gJTT9G+HV/XtqpqvqkvdxScJRX30vzV7uOn1+bSOb8Antw+xhMIYY4xXzu/enA8mDOZYsXLFS3NYsOVgoEMqkz3EryHTlu9i/BsL6dQsmnfHDzqtq6+MMcacfrolxfLxr4YQHxXGda/+xNcrdwc6pJ+xpKIGfLZkBxPfXkSvVnG8detAGyHTGGNMlbSKb8CHEwbTqXkMv3prIe/P3175STWo0qTC7XnxCxF5wF1vLSID/B9a3fD1yt3c8/5SBrSNZ/LNA4iJCA10SMYYY2qxhIbhvPPLgQxNSeQPHy3j7Z+2BTqk47ypqXgBGIwz6BVAFs68G6YSM9fv4863F9OjZSz/GdefqHBvOtsYY4wxFWsQFsLL1/dlRGoi93+ynDfnbg10SIB3ScVAVb0DyANwh7y2+vtKLNx6kPFvLKRdYhSv3zjAEgpjjDE+FREazEvX9+Xszk3486crmDx7S6BD8iqpOObOCqoAIpIIFPs1qlpuy/4cbp28gGaxEbx5y0BiG9gjD2OMMb4XHhLMC9f1ZXSXpjw4dSUfL8oIaDzeJBXPAp8ATUTkMWAW8H9+jaoWO5J7jJtfnw/Aazf2JzE6PMARGWOMqcvCQoJ47trenNEhgd9/uIzvVu+p/CQ/qTSpUNUpwB+AvwK7gEtU9QN/B1YbFRQWM+GthWw/lMu/r+9HcuOoQIdkjDGmHggPCebf1/eja4sYbp+yiHmbAzOOhbddStfj1FZMBXJEpLX/Qqq9/m/aauZsOsATl/dgQNv4QIdjjDGmHmkYHsJrN/YnqVEkv3xjAVv259R4DN50Kb0T2AN8izPj5xecmPnTuL5YtovXZ2/h5jPaclmfloEOxxhjTD2U0DCc128cQJDALZPnk5l3rEbv701NxV1Aqqp2VdUeqtpdVXv4O7DaZNO+bP740TJ6t47j3vM6VX6CMcYY4yetExrwwnV92XoglzvfXkxRsdbYvb1JKrYDR/wdSG1VUFjMxLcXExosPH9tH8JCbJBSY4wxgTW4fQKPjOnGD+v28dTXa2vsvuUOniAi97hvN+FMff4FcHwyMVV92s+x1Qr/+t96Vu3K5N/X96VFXGSgwzHGGGMAuHZga5bvOMxLP2xkcPsEzuyY6Pd7VvS1OtpdtuG0pwjz2NbQ75HVAssyDvN8+kYu653EOV2bBTocY4wx5iQPXtSV1KbR3PPeEvZk5vn9fuUmFar6sKo+DKwqee+xbbXfIzvNHSsq5ncfLKVxwzAevKhroMMxxhhjfiYiNJh/XdubnIJC7n5vCcV+bl/hTQOA+7zcdspE5FwRWSsiG0Tk3jL2i4g86+5fJiJ9KjtXRJ4SkTXu8Z+ISJy7PVlEjorIEnd5qTqxT569hXV7svnLmG42YqYxxpjTVkrTaB68qCuzNx5gyjz/Tj5WUZuK84DzgSQRedZjVwxQWN0bu0N/Pw+MAjKA+SIyVVVXeRx2HpDiLgOBF4GBlZz7LXCfqhaKyBM4CdAf3ettVNVe1Y19b2Yez0xfT1pqIqO6NK3u5Ywxxhi/uqZ/K75YtovHp63mrE5NSPJTG8CKaip2AgtwJhJb6LFMBc7xwb0HABtUdZOqFgDvAmNKHTMGeEMdc4E4EWle0bmq+o2qliQ9cwGfDxrxxFdrKSgs5sGLuiIivr68McYY41Miwl8v644C9328HFX/PAYpt6ZCVZcCS0XkbVX1x+gZSTjdVUtk4NRGVHZMkpfnAtwMvOex3lZEFgOZwP9T1ZllBSYi44HxAImJiaSnpx/ftyO7mI8XHeWc5FC2rpjP6THZbO2WnZ19Uhkb37MyrhlWzv5nZVw9l7YPZsrqffz9ve/o18z3s2dXekU/JRQAZX3FL506lXdMpeeKyJ9wHtNMcTftAlqr6gER6Qt8KiJdVTXzZxdSfRl4GSA1NVXT0tKO77tjyiIahBXwf9enER9lM8D7Qnp6Op5lbHzPyrhmWDn7n5Vx9QwtKmbBs7P4bFshEy8fRkRosE+vH8iRmjKAVh7rLXEeuXhzTIXnisg44ELgOnXreFQ1X1UPuO8XAhuBjqcS8KqdmXyxfBe3DG1rCYUxxphaJyQ4iD9f2IXtB48y6cfNPr9+uUmFiLzpvt7l87s65gMpItJWRMKAa3Daa3iaCtzg9gIZBBxR1V0VnSsi5+I0zLxYVXM9Pk+i28ATEWmH0/hz06kE/MrMTUSFBXPLsHZV+bzGGGNMwA1NaczZnZvywvcbOZLr24cRFdVU9BWRNsDNItJIROI9l+re2G1MORH4Gmfci/dVdaWITBCRCe5h03D+8G8AXgFur+hc95x/4QzQ9W2prqPDgWUishT4EJigql7PDbsnM4//LtvJlf1aERtpXUiNMcbUXr8d3ZHs/EL+4+PaioraVLwEfAW0w+n14dmOQd3t1aKq03ASB89tL3m8V+AOb891t3co5/iPgI+qGutbc7dSWKzcdEZyVS9hjDHGnBY6N4/hnK5Nee3HzdwytK3PvixXNKLms6raGZikqu1Uta3HUq/q/4uKlfcXbGdEahPaJEQFOhxjjDGm2n49MoWsvELemuu7foyVNtRU1V+JSE8Rmegu9W7a8zkbD7AnM5/L+/h8yAtjjDEmILq2iGVI+wTe/mmbz6ZHrzSpEJFf43TLbOIuU0TkTp/cvZb4ZPEOosNDGNm5SaBDMcYYY3zm+kFt2HH4KN+t3uOT63kz8sWtwEBVzQFwh76eAzznkwhOcwp8s2o353Rr5vP+vMYYUx3Hjh0jIyODvDz/zz55uoiNjWX16no/p6VPREREkJbSgmYxEbw5dyujfTDbtjdJhQBFHutFlD34VJ2UXwQFeYWcbbUUxpjTTEZGBtHR0SQnJ9ebKQOysrKIjo4OdBi1nqpy4MABdu/ayeV9k3gxfSP7s/Np3DC8Wtf1ZvCr14CfROQhEXkIZz6N/1TrrrXI0WNKSJBwRofGgQ7FGGNOkpeXR0JCQr1JKIzviAgJCQnk5eVxUc8WFCt8uWJ3ta/rTUPNp4GbgIPAIeAmVX2m2neuJY4WKn1aNyI6wsamMMacfiyhMFVV8rOT2jSalCYN+e/S0oNanzqvhulW1UVuF9N/quriat+1Fikohr7JjQIdhjHGnHa2bNlCt27dTtr20EMP8be//S1AEfnXAw88wPTp031yrbFjx9KjRw/+8Y9/+OR61SEijO7alAVbD5GVV70RNn0/RVkd1LNlXKBDMMYYg9MWoLi4mKAg309dpaqoarnXfuSRR3xyn927dzN79my2bq3a+BCFhYWEhISUu+7teZ7O6NCY57/fyNxNXg80XaZATihWa/RsFRvoEIwxptZJS0vjN7/5DUOGDKFbt27MmzcPcGozrr/+es466yxSUlJ45ZVXjp/z1FNP0b9/f3r06MGDDz4IODUinTt35vbbb2fYsGFs3779pPskJydz//33M3jwYPr168eiRYs455xzaN++PS+95AzSnJ2dzciRI+nTpw/du3fns88++9m1+/Tpw/bt2/nLX/5Cp06dGDVqFGPHjj1e83LjjTfy4YcfHr/ngw8+ePx6a9as+dnnz8vL46abbqJ79+707t2b77//HoDRo0ezd+9eevXqxcyZM086Z9++fVx++eX079+f/v378+OPPx4vs/HjxzN69GhuuOGGn61v3bqVkSNH0qNHD0aOHMm2bduOx3zPPfcwYsQI/vjHP5b7b9W3TSMiQoOYtX6fN/+05aowtXEn4PpaVc+u1l1qMQGaxUQEOgxjjKnQw5+vZNXOTJ9es0uLGB68qGu1rpGTk8Ps2bOZMWMGN998MytWrABg2bJlzJ07l5ycHHr37s0FF1zAihUrWL9+PfPmzUNVufjii5kxYwatW7dm7dq1vPbaazzxxBNl9v5o1aoVc+bM4e677+bGG2/kxx9/JC8vj65duzJhwgQiIiL45JNPiImJYf/+/QwaNIiLL74Y4Pi1X3jhBRYsWMBHH33E4sWLKSwspE+fPvTt27fMz9a4cWMWLVrECy+8wN/+9jdeffXVk/Y///zzACxfvpw1a9YwevRo1q1bx9SpU7nwwgtZsmTJz6551113cffddzN06FC2bdvGOeecc7wL7cKFC5k1axaRkZE89NBDJ61fdNFF3HDDDYwbN45Jkybx61//mk8//RSAdevWMX36dIKDyx8WITwkmN6tGrF4++EK/z0rU2FSoapFIpIrIrGqeqRad6qlQoOsIZQxxpSlvN+NntvHjh0LwPDhw8nMzOTw4cMAjBkzhsjISCIjIxkxYgTz5s1j1qxZfPPNN/Tu3RtwahfWr19P69atadOmDYMGDSIrK6vMe5YkCN27dyc7O5vo6Giio6OJiIjg8OHDREVFcf/99zNjxgyCgoLYsWMHe/Y4Az6VXBtg1qxZx2MDuOiii8r9/JdddhkAffv25eOPP/7Z/lmzZnHnnc5YkZ06daJNmzasW7eOmJiYcq85ffp0Vq1adXw9MzPz+Ge++OKLj8dVen3OnDnHY7j++uv5wx/+cPy4K6+8ssKEokT3lrG8PntLpcdVxJs2FXnAchH5Fsgp2aiqv67WnWuJkCBLKIwxp7/q1ihURUJCAocOHTpp28GDB2nbtu3x9dKJR8l6WdtVlfvuu4/bbrvtpH1btmwhKqrieZfCw53xFYKCgo6/L1kvLCxkypQp7Nu3j4ULFxIaGkpycvLxQcM8r+3MY+mdkvsEBwdTWFj4s/2ncq0SxcXFzJkz56TkoUTpMqioTDzLt7KyK9G1RQwFhcVeRlo2b9pUfAH8GZiBM1tpyVIvBFtOYYwxZWrYsCHNmzfnu+++A5yE4quvvmLo0KHHj3nvvfcA51t7bGwssbFOG7XPPvuMvLw8Dhw4QHp6Ov379+ecc85h0qRJZGdnA7Bjxw727t3rk1iPHDlCkyZNCA0N5fvvvy+3keTQoUP5/PPPycvLIzs7my+++KLK9xw+fDhTpkwBnEcQ27ZtIzU1tcJzRo8ezb/+9a/j62U9IinLkCFDePfddwGYMmXKSf8G3urcvPwaFG9VWlOhqpNFJBJoraprq33HWsYqKowxpnxvvPEGd9xxB7/97W8BePDBB2nfvv3x/Y0aNWLIkCFkZmYyadKk49sHDBjABRdcwLZt2/jzn/9MixYtaNGiBatXr2bw4MGAk7S89dZbXlXdV+a6667joosuol+/fvTq1YtOnTqVeVz//v25+OKL6dmzJ23atKFfv37HE6FTdfvttzNhwgS6d+9OSEgIr7/++km1KGV59tlnueOOO+jRoweFhYUMHz78eGPTys67+eabeeqpp0hMTOS111475XhbNvp57cgpK+lCU94CXASsBTa7672AqZWdV1eWJm1S1PjX999/H+gQ6jwr45pR0+W8atWqGr3fqTrzzDN1/vz5P9v+4IMP6lNPPVWla2ZmZlY3rEplZWWpqmpOTo727dtXFy5c6Pd7Bkrpn6Hej3yjwAKt4t9Mbx5/PAQMAA67ScgSoG35h3tPRM4VkbUiskFE7i1jv4jIs+7+ZSLSp7JzRSReRL4VkfXuayOPffe5x68VkXO8iTHEOt0aY0y9Mn78eHr16kWfPn24/PLL6dOnT+Un1RFNq9nb0ZuGmoWqeqRUo5pqT7zudld9HhgFZADzRWSqqq7yOOw8IMVdBgIvAgMrOfde4DtVfdxNNu4F/igiXYBrgK5AC2C6iHRUVc/J0n6mYag9/zDGmKpIT08vc/tDDz1Uo3GcqrfffjvQIQRMTET1xsT05nv4ChG5FggWkRQReQ6YXa27OgYAG1R1k6oWAO8CY0odMwZ4w62VmQvEiUjzSs4dA0x2308GLvHY/q6q5qvqZmCDe50KhRTmVHaIMcYYUyfERFZvnitvUpI7gT8B+cDbwNfAX6p1V0cS4DksWgZObURlxyRVcm5TVd0FoKq7RKRkzvIknBlWS1/rZ0RkPDAeoEPThuVm28Y3srOzrYz9zMq4ZtR0OcfGxpY7bkNdVVRUVO8+sz/l5eWd9DMbcrSgWtfzJqm4QFX/hJNYACAiVwIfVOvOzmCVpZV+rFLeMd6cW5X7ORtVXwZeBuiRnKBpaWmVXNpUR3p6OlbG/mVlXDNqupxXr15d5uiSdVlWVla9+8z+FBERcXywMYC0NHhpQtWv583jj/u83HaqMoBWHustgdLzrpZ3TEXn7nEfkeC+lnRy9uZ+PyNavYFAjDHGmFpjzgvVOr3cpEJEznPbTyS5PTBKlteBnw8ddurmAyki0lZEwnAaUU4tdcxU4Aa3F8gg4Ij7aKOic6cC49z344DPPLZfIyLhItIWp/HnvMqClIrbcRpjTL3WsGHDQIdQY4YMGeKT66xZs4ZevXrRu3dvNm7c6JNr+sy2OdU6vaLHHzuBBcDFnDyCZhZwd7XuCqhqoYhMxGmjEQxMUtWVIjLB3f8SMA04H6dRZS5wU0Xnupd+HHhfRG4BtgFXuuesFJH3gVU4SdEdlfX8AAgq9kX+ZIwxxhfKGg7bl9euaArx2bN90UcBPv30U8aMGcPDDz9cpfOLiopOGhCs9Lq355Xp2NEqxVSi3JoKVV2qqpOBDsD7wFxVnayqH6vqofLOOxWqOk1VO6pqe1V9zN32kptQ4Pb6uMPd311VF1R0rrv9gKqOVNUU9/Wgx77H3ONTVfVLb2IM0mNQbI9AjDHGW0uWLGHQoEH06NGDSy+9lEOHDrF3797js30uXboUETk+PXf79u3Jzc31atrv8ePHn3Sv9PR0zjzzTK666io6duzIvffey5QpUxgwYADdu3c/XhPw+eefM3DgQHr37s3ZZ599fDKx0lOI79u3j1GjRtGnTx9uu+022rRpw/79+4ETtTIlbWeuuOIKOnXqxHXXXVfmPB9llcO0adN45plnePXVVxkxYsTPzvnmm28YPHgwffr04corrzw+ZHlycjKPPPIIQ4cO5YMPPvjZ+jvvvEP37t3p1q3bSVOcN2zYkAceeICBAwcyZ44XtRD51Zvp1puGmucCfwPCgLYi0gt4RFUvrtadawsthqxdEFtmRxFjjDk9fHkv7F7u22s26w7nPX7Kp91www0899xznHnmmTzwwAM8/PDDPPPMM+Tl5ZGZmcnMmTPp168fM2fOZOjQoTRp0oQGDRpw6623Vjrtd1k1FUuXLmX16tXEx8fTrl07br31VubNm8c///lPnnvuOZ555hmGDh3K3LlzERFeffVVnnzySf7+97+fdO3IyEgmTpzIWWedxX333cdXX33Fyy+/XOZnXLx4MStXrqRFixacccYZ/Pjjjz+bb6O8cpgwYQINGzbkd7/73UnH79+/n0cffZTp06cTFRXFE088wdNPP80DDzwAOI0qZ82aBcC99957fH3nzp0MGjSIhQsX0qhRI0aPHs2nn37KJZdcQk5ODt26deORRx7x7h/vyA7vjiuHN0nFQzjjOaSDM6KmiCRX6661zZ4VllQYY4wXjhw5wuHDhznzzDMBGDduHFdeeSXgtEn48ccfmTFjBvfffz9fffUVqsqwYcMA76b9Lqs7af/+/WnevDng1HqMHj0acKZB//777wHIyMjg6quvZteuXRQUFJw0k6rnFOKzZs3ik08+AeDcc8+lUaNGlGXAgAG0bNkSgF69erFly5aTkoqKyqE8c+fOZdWqVZxxxhkAFBQUHJ8HBeDqq68+6fiS9fnz55OWlkZiYiLgzHMyY8YMLrnkEoKDg7n88ssrvO9xxUXOl+hqqOqImvWIwPZ50NGrUb2NMSYwqlCjUNOGDRvGzJkz2bp1K2PGjOGJJ55ARLjwwguBU5v221Ppqc49p0Evqdm48847ueeee7j44otJT08/aVTPqkx97nnP8qY+P1WqyqhRo3jnnXfK3F/e1OcVxRwREeH9hGyZO6GanRMCOaJmrVAcHAbbfwp0GMYYUyvExsbSqFEjZs6cCcCbb755/Nv68OHDeeutt0hJSSEoKIj4+HimTZt2/Jt5Vaf99saRI0dISnJqnCdPnlzucUOHDuX9998HnPYNhw5VrQlhReVQnkGDBvHjjz+yYcMGAHJzc1m3bl2l9xo4cCA//PAD+/fvp6ioiHfeeafSe5Vpz8rKj6nEqY6o+Q6+G1GzVigMjoSMBU6L2FAfTAtrjDF1SG5u7vHHAAD33HMPkydPZsKECeTm5tKuXbvj03AnJycDTnIBzh/wjIyM448YqjrttzceeughrrzySpKSkhg0aBCbN28u87gHH3yQsWPH8t5773HmmWfSvHnzKg+2VV45lCcxMZHXX3+dsWPHkp+fD8Cjjz5Kx44dKzyvefPm/PWvf2XEiBGoKueffz5jxpSe9cILe6rfJke8reqpr7q2b6Urr8+E6z6ElFGBDqdOstEe/c/KuGYEYkTNzp0719j9Tgf+HlEzPz+f4OBgQkJCmDNnDr/61a98WmNyujnpZ+jd62DPSuQ3Sxeqar+qXK/SmgoR6QfcDyR7Hq+qPapyw9qmKCQSQgth3VeWVBhjTB23bds2rrrqKoqLiwkLC+OVV14JdEg1o7gYtsyCzhcCS6t8GW8ef0wBfg8sB+rdgA2KQPuzYM00OO9JCPKywYsxxphaJyUlhcWLFwc6jJq3ZznkHYbk4UDVh+r2pqHmPlWdqqqbVXVryVLlO9ZG3a+ArJ2w+YdAR2KMMcb43sb/Oa9th1XrMt7UVDwoIq8C3+E01gRAVT+u1p1rk47nQUQsLHnHqbUwxpjThKpSf7v8m+o4qU3lqqnQojfEtKjWNb1JKm4COgGhnHj8oUD9SSpCI6DbFbDkbcg9CA3iAx2RMcYQERHBgQMHSEhIsMTCnBJV5cCBA0RERMChrbBzEZz9ULWv601S0VNVu1f7TrVd/1tgwX9g4Wsw7LeBjsYYY2jZsiUZGRns27cv0KHUmLy8POcPoam2iIgIpzvwT887G7pUoRtqKd4kFXNFpIuqrqr80DqsaVdoNwJ+ehkG3wkhYYGOyBhTz4WGhp403HR9kJ6eTu/evQMdRt1RXAyLJkPrwRDfrtqX86ah5lBgiYisFZFlIrJcRJZV+8610ZCJkL0blr8f6EiMMcaY6tv8AxzcBP1u9snlvJ2l1AC0HwnNe8IPT0D3q6y2whhjTO02/1WIjIfOvpl4vNKaCo8upEdxGmiWLPWPCIx8AA5vc6qLjDHGmNpq31pY8wX0u8npkOADlSYVInKxiKwHNgM/AFuAL6tzUxGJF5FvRWS9+1rm3LIicq772GWDiNxb2fkiMkpEFrqPaBaKyFke56S711riLk2qFHz7kdDmDPjhScjPrtIljDHGmICb+XdnTqtBd/jskt60qfgLMAhYp6ptgZHAj9W8773Ad6qagjP+xb2lDxCRYOB54DygCzBWRLpUcv5+4CK3t8o44M1Sl71OVXu5y94qRS4Cox6BnL3OYxBjjDGmtjmwEZZ/4LSliErw2WW9SSqOqeoBIEhEglT1e6BXNe87Bih5fjAZuKSMYwYAG1R1k6oWAO+655V7vqouVtWd7vaVQISInJj03lda9oPe18PcF2Dvap9f3hhjjPGr6Q9BSAQMudOnl/WmoeZhEWkIzACmiMheoLCa922qqrsAVHVXOY8ikoDtHusZwMBTOP9yYLGq5ntse01EioCPgEe1nClaRWQ8MB6cqWjT09N/dkxo5CgGBH1CzpRbWdLrUacGw1RJdnZ2mWVsfMfKuGZYOfuflXH1xR5eSe/VU9mcfC1bF64B1vjs2t4kFWNwGmneDVwHxAKPVHaSiEwHmpWx609exlbWX2mvGoiKSFfgCWC0x+brVHWHiETjJBXXA2+Udb6qvgy8DJCamqrlTmWccIi4z+8iLWoDDPilN6GZMti03P5nZVwzrJz9z8q4moqL4dWHILoFba/9O23DGvj08hUmFW67hs9U9WycIbq97vLgnlPedfeISHO3lqE5UFb7hgyglcd6S6Dk0Ua554tIS+AT4AZV3egRzw73NUtE3sZ5vFJmUuG1PuOc8dK/+bMzJ0hC+2pdzhhjjPGr+a/CzsVw6cvg44QCKmlToapFQK6IxPr4vlNxGlLivn5WxjHzgRQRaSsiYcA17nnlni8iccAXwH2qerwxqYiEiEhj930ocCGwotqfQgTG/MsZr+KT26Couk+FjDHGGD85vB2+e9j5EtzjKr/cwpuGmnnAchH5j4g8W7JU876PA6Pcrqqj3HVEpIWITANQ1UJgIvA1sBp4X1VXVnS+e3wH4M+luo6GA1+7I4EuAXYAr1TzMzhiWsAFT0PGfPj+MZ9c0hhjjPEpVfjiHuf1wmf81g7QmzYVX7iLz7i9SUaWsX0ncL7H+jRg2imc/yjwaDm37VvVeCvV/QpnqNNZT0OrAZB6nt9uZYwxxpyyha/B+m/g3MehURu/3abSpEJVbehIb5z3FOxaCh/fBrf9APH1a5IfY4wxp6k9q+Cr+5zBGwfc5tdbeTOiZoqIfCgiq0RkU8ni16hqo9AIuOoNp8/Ku9dCXmagIzLGGFPfFeTChzdBeAxc+hIEedPqoeq8ufprwIs4Y1OMwOkxUXqkSgPQKBmunAz718EHN1rDTWOMMYGjCv+9G/atgcv+DQ2rNjvFqfAmqYhU1e8AcScXewg4q5Jz6q/2I5yGmxu/gy9/7/yjGmOMMTVtzr9g2buQdr/T46MGeNNQM09EgoD1IjIRp+eE/9Od2qzvOGd++h+fgdiWMOy3gY7IGGNMfbL+W/j2AehyCZz5hxq7rTdJxW+ABsCvcSYXO4sTY0SY8ox8EDJ3wnePQFg0DBwf6IiMMcbUB3tWwoe3QNOucMkLNTqNhDe9P+YDuLUVv1bVLL9HVRcEBTn/mAU5zmOQsCjofV2gozLGGFOXHdoKb17m/M255h3ntQZ50/ujn4gsB5bhDIK1VET8N+ZDXRIcCle+Bu1GwNSJsPzDQEdkjDGmrsreB29eCoV5cP3HENeq8nN8zJuGmpOA21U1WVWTgTtweoQYb4SEwzVToPUQ+OhWWPxWoCMyxhhT1xw9DFOucB67X/s+NOkckDC8SSqyVHVmyYqqzgLsEcipCIuC6z5weoZ8dgfM880I4cYYYwxHD8GblzhtKa56A1oPDFgo3iQV80Tk3yKSJiJnisgLQLqI9BGRPv4OsM4IawBj34XUC2Da72DWP6y7qTHGmOrJPQhvjHESiqvfgo6jAxqON70/ermvD5baPgRQbMwK74WEw1WT4ZMJMP0hp5rq3MchKDjQkRljjKltShKKfWvh6ikBTyjAu94fI2oikHojOBQuewVimsPs5+BIBlz+ao230DXGGFOLHdoKb10Oh7fBNW9DytmBjgjw7vGH8bWgIBj9KJz/N1j3Fbx+IWTvDXRUxhhjaoPdy+E/oyFnL9zw2WmTUIAlFYE14JdOldXe1fByGuxYFOiIjDHGnM42z4TXzncem9/8NbQZHOiITmJJRaB1Oh9u+RokGCadC0veDnRExhhjTkcLJzvjUEQ3h1u+CVi30Yp4M/hVAxH5s4i84q6niMiF1bmpiMSLyLcist59bVTOceeKyFoR2SAi91Z2vogki8hREVniLi95nNNXRJa713pWpAbHLa1M854wPt3pBvTpr2DaH6DoWKCjMsYYczooKnT+Lnz+a2g7zPkiGtsy0FGVydupz/OBkjqWDODRat73XuA7VU0BvnPXTyIiwcDzwHlAF2CsiHTx4vyNqtrLXSZ4bH8RGA+kuMu51fwMvhWVAL/4BAZPhHn/htcvcBrgGGOMqb9yD8Jblzl/FwZPhGs/gMgyv4efFrxJKtqr6pPAMQBVPQpU91v+GGCy+34ycEkZxwwANqjqJlUtAN51z/P2/ONEpDkQo6pzVFWBNyo7JyCCQ+Ccx+CKSbBnFbw0FFZNDXRUxhhjAmHHIqe93bY5MOYF5+9DsDcjQQSON9EViEgkzpgUiEh7nJqL6miqqrsAVHWXiJQ1lXoSsN1jPQMoGSasovPbishiIBP4f+5ooEnu+Z7XSiovOBEZj1OrQWJiIunp6afy2XwggYjeT9Fl1d+Jef96drQ4j43tb6I4OLyG46gZ2dnZASjj+sXKuGZYOftfvShjVZJ2fEH7ja9REBbHqh6PknkkCWrB5/YmqXgQ+ApoJSJTgDOAGys7SUSmA83K2PUnL2MrqzaksiEodwGtVfWAO+nZpyLS9VSvpaovAy8DpKamalpamncR+9qoy+F/j5A0+zmSCrfCpS857S/qmPT0dAJWxvWElXHNsHL2vzpfxnlH4LOJsGEqdDyXiEtepE+D+EBH5TVvBr/6VkQWAYNw/jjfpar7vTiv3I6zIrJHRJq7tQzNgbIGacgAPKdYawnsdN+Xeb6q5uPWoqjqQhHZCHR0r9WynGudvkLCnPEs2p7pzBnyylkw/Pcw7LfOIFrGGGPqjm1z4ePxkLnD+d0/eCKcRn0KvOFtl9II4BDOI4UuIjK8mvedCoxz348DPivjmPlAioi0FZEw4Br3vHLPF5FEt4EnItIOp0HmJvdRSZaIDHJ7fdxQzj1PTymj4Pa50PVSSP8rvDrSaXNhjDGm9ivMh28fcIYVQOHGaTDkzlqXUIAXNRUi8gRwNbASKHY3KzCjGvd9HHhfRG4BtgFXuvdqAbyqqueraqGITAS+BoKBSaq6sqLzgeHAIyJSCBQBE1T1oLvvV8DrQCTwpbvUHg3ineG8O18E/70HXj4TzvwDDLnLqdEwxhhT++xa5swHtXcl9BnnNMYMjw50VFXmTZuKS4BU99GCT6jqAWBkGdt3Aud7rE8Dpp3C+R8BH5VzzwVAt6pHfZroMgbanAFf/Bb+9ygs/xAu/Ae0GRLoyIwxxnirsABm/xPSn3C6iF77PnQ8J9BRVZs3jz82AfYA/3QS1diZ7XTse1CQC6+d57S5yD1Y+bnGGGMCa/s8+Pdw54thpwucx9t1IKGACmoqROQ5nMccucASEfkOj66kqvpr/4dnKpR6rjO62g9PwOx/wdovYdQj0PNaZ9IyY4wxp4+8I/DdIzD/PxCTBGPfhdTzAh2VT1X0+GOB+7qQEw0kS1TWtdPUlLAoJ5HocTV8/hunxmL+q3Du49B6UKCjM8YYowqrPoOv7oXsPTBwApz1p1rddqI85SYVqjoZQETuUtV/eu4Tkbv8HZg5RU27OjPWLf8Apj8Ek86BbpfD2Q9DXKtKTzfGGOMHu1c4ycSWmdCsO1wzBZL6Bjoqv/GmjnxcGdtu9HEcxheCgqDn1XDnAjjzj7DmC/hXP/jfY5CfFejojDGm/sg5AP+9G/49DPashAuehl+m1+mEAipuUzEWuBZn2GvPxx/RwAF/B2aqISwKRtwPva93ai1mPAkLJjmDZvW7GUIjAh2hMcbUTYUFsOA/zphC+dkwYLzzJa8WjYpZHRW1qZiNM+x1Y+DvHtuzgGX+DMr4SFwruOI/MOh2+O5h+Po+mPM8pN0LPcee9hPTGGNMrVFcDCs+dHp0HN4K7dKctm1NOgc6shpVUZuKrcBWTkx5bmqrln1h3FTYlA7TH4apE2H2szDiT9D5YuspYowxVaUKG6Y7v1v3LHfaTVz3EXQYWStHxKwu+2tSn7RLg1/+D656ExD4YBy8dIYzgFZxUaCjM8aY2mXbT/D6BTDlCijIgsv/A+NnQMrZ9TKhAEsq6h8R6HIx/Go2XPaKk0x8dAs8PwCWvA1FhYGO0BhjTm9bZsHki2HSaDiwAS74O9wxH7pfUe9rfuv3p6/PgkOgx1XOSG5XToaQSPj0V/BcH1j4OhzLC3SExhhz+lCFjd/Da+c7tRP71sDox+DXi6H/rTYHk8ubCcXOAB4C2rjHC6Cq2s6/oZkaERQEXS9x5hRZ+6XTU+Tzu5xuqAPGO71FohICHaUxxgSGKqz/FmY8BRnzILoFnPck9LkBQiMDHd1px5vm//8B7sYZWdMevNdVItDpfGfI2E3pMOdf8P2jMPPv0GssDLoDGncIdJTGGFMzjuXB8vedHnP71kBsK2esid6/gJDwQEd32vImqTiiqrVrmnBTdSLQfoSz7F3t/Ida/BYseM1JOAbdDslD620jJGNMHZd70JmbY97LkLMXmnaHS/8NXS+zRxxe8Cap+F5EngI+5uQJxRb5LSpzemjSGcb8C0Y+APNeceYUWTsNEjs5zxB7XA0RMYGO0hhjqm/fOpj3b1g8BQqPQoezYcid0PZM+xJ1CrxJKga6r/08tilwlu/DMaelhk2cyW+G3QMrPnaSi2m/g28fdBp79r8VmnULdJTGGHNqio45X5TmvwqbZ0BwmPM7bfDEejdola9UmlSo6ghf31RE4oH3gGRgC3CVqh4q47hzgX8CwcCrqvp4ReeLyHXA7z0u0QPoo6pLRCQdaA4cdfeNVtW9vv5sdVpoJPS+zll2LIT5k2DpO7DwNWg1CPrd5AymFdYg0JEaY0z5MnfBoslOT7esXU57iZEPQO8boGFioKOr1SrtUioisSLytIgscJe/i0hsNe97L/CdqqYA37nrpe8bDDwPnAd0AcaKSJeKzlfVKaraS1V7AdcDW1R1icdlryvZbwlFNSX1hUueh3tWO92qcvbCJ7fB3zrC1F/D9nlOq2ljjDkdFBfB+unw3vXwj67O3BxNu8LYd+Gupc7cSJZQVJs3jz8mASuAq9z164HXgMuqcd8xQJr7fjKQDvyx1DEDgA2quglARN51z1vl5fljgXeqEaPxRoN4GDLRacC5bbbzPHL5B863gMYdode1zjwj0c0CHakxpj7avwGWTHFqVbN2QWQ8DL7d6S4fbyMj+JpoJd8mRWSJ+82/wm2ndFORw6oa57F+SFUblTrmCuBcVb3VXb8eGKiqE708fyMwRlVXuOvpQAJOt9iPgEe1nA8vIuOB8QCJiYl933///ap+1HopuDCXxH0/0nzXd8RmrkYJ4mB8b/Y0Hc6BhIEUhZzctzs7O5uGDRsGKNr6wcq4Zlg5+583ZVz276A+7Go+kgMJ/dGg0BqKtnYaMWLEQlXtV/mRP+dNTcVRERmqqrPg+GBYRys5BxGZDpT19fRPXsZWVnNbr+rTRWQgkFuSULiuU9UdIhKNk1RcD7xR1vmq+jLwMkBqaqqmpaV5GbI54XzgMdi/AVnyFgnLPiBh9T+ckTtTz4XuVzqtq0PCSU9Px8rYv6yMa4aVs/+VW8aFBbDxO2cuo7XT4FiuU1t69sNIz2tIiG6GDePnf94kFROANzzaURwCxlV2kqqeXd4+EdkjIs1VdZeINAfKat+QAbTyWG8J7HTfV3b+NZR69KGqO9zXLBF5G+fxSplJhfGhxh3g7IfgrAdg+1znP/yqT2HlJxARC50vIq64IxQPg6DgQEdrjKlNiotgy0zn98rqqZB3BCIbOd3de10HLftZd9Aa5k3vj6VATxGJcdczfXDfqTiJyePu62dlHDMfSBGRtsAOnETh2srOF5Eg4EpguMe2ECBOVfeLSChwITDdB5/DeCsoCNoMcZbznoBNPzhtL1Z+Sq+CbFj/rDOiZ6eLoN2ZNmKdMaZsWuzMDrryY+fLSfYeCGsInS6Ablc4A/cF2+ONQPGmpgLwWTJR4nHgfRG5BdiGkwQgIi1wuo6er6qFIjIR+BqnS+kkVV1Z0fmu4UBGSQNPVzjwtZtQBOMkFK/48POYUxEc6kwNnHI2HDvKyk+epmvQJljxCSx6A8KioeM50PlC6DAKwu0ZtTH1WmGBUyOx5r8MXvYJ/HAIgsOh42jodjmknGNd2U8TXicVvqSqB4CRZWzfifMwvmR9GjDN2/PdfenAoFLbcoC+1Qra+EdoJPuanAFpf4LCfKcGY83nsOYLWPGh84uj/VnOEOEpoyGmeaAjNsbUhIIc2DAdVv8X1n0N+UcgtAGZsT1JHHaT0zYrorqjGxhfqzCpcB8lDFLV2TUUj6nPQtxvHh1HwwX/cNpgrP7c/aXiTj/TrLvzrSRltPO81NphGFN3HN7mzAi6/lvY9D0U5jltJDpfCJ0uhPYjWPnjT6T1TAt0pKYcFSYVqlosIn8HBtdQPMY4gkOcicuSh8K5j8PeVc63lfXfwqx/wMy/Ob9sOpztJBjtR9oU7cbUNoUFsG0OrP/GqZXYt8bZHtca+oxzkonWQ5zfB6ZW8OZf6hsRuRz4uLxxHYzxKxFn5LumXZ35R44ego3/O/GNZvkHgDi1GO3SnIaerYfYM1ZjTjeqcGgLbEp3kohN6VCQ7cy50WYI9LnBaUfVOMV6bdRS3iQV9wBRQJGIHMUZP0JV1aanNIER2chpnNXtciguhp2LYOP3zi+ouS/C7GedX1ItB7hJRhq06G3fdowJhKzdzmRdm3+ATTPgyDZne2wrZ/KuDqOg7XBrkF1HeNOlNLomAjGmSoKCnLYVLfvBmb93Gndtm+M0+NyUDt8/Bt8/CuEx0GogtBns1GIk9bFuq8b4Q+5B2Pqj839w8wzYv9bZHhEHbYfBGb92phO32og6yauvbiJyMSfGfUhX1f/6LyRjqiEsymln0cEdey3nAGyZ4fyC2zYHvvvW2R4c7kyK1nqQU+3aaoC1JDfmVJU8ztg21/n/tW3uiSQitIHzf6v3L5yaiGbdrWF1PVBpUiEijwP9gSnuprvcYbt/NrOoMaedqAToeqmzgJNkbJ8LW2c7vwRnPwuzngYJctpstOwPSf2chKNxR6cmxBjjKDoGu5fD9p9OJBHZe5x94bHQeqDzSKPNGc7/oZCwwMZrapw3NRXnA71UtRhARCYDiyljunJjTntRCc7Ie50ucNYLciBjPmyd4/ySXP4hLJjk7AuPcdpiJPV1Hq8k9bXZVk39UVwMBzfCjkVOu6Udi2D3MqebJzg9NNqe6dT2tR4MiZ0sCTdeD34VBxx031sdsak7wqJONOYE5xfpgfWQsQB2LIQdC5zajOJCZ39MS2jRC5r1cKpzm3WH2Jb2bNjUbqpwZDvsXHIigdi5xBlwCpxHGc17Qr9boGVfaDUIYpMCGbE5TXmTVPwfsFhEvsfp+TEcuM+vURkTKEFBkJjqLL2vc7YdOwq7lp1IMnYtc0b8LJk0N7KRm2B4JBqNO9r8A+b0dOwo7F0Ne1bA7hXO654VzmRcAEEhzqPAbpc5tXNJfaBxqvWeMl7xZkTNYpxhr/vjJBV/VNXdNRCbMaeH0EjnWXHrgSe2FeTAnpVOdfCuZc5z5vmvnqgaDg6DhBQ3QekETTo5r/HtLNkwNaOo0GlEuX+tM6jUnpVOEnFgvTMpF0BoFDTtAl0vg2bdoHkvaNoNQiMCGbmpxbwZUXOiqr6PMzOoMQacxyatBjhLiaJC5xf27uXOsm+tU7ux8uMTxwSFQkIHJ9lo0tnpVhff3kk2ImzoF1MFBbnOz92+dU4CsX+d8/7gRigqOHFcbGsncegyxqmJaNYdGrW1dhDGp7ypz/pWRH4HvAfklGxU1YPln2JMPRQc4iQKTTo7LeBLFOQ4v+j3rnG+Me5bC7uWwKrPOP4IBSAq0Uku4ttDQrsT7y3hMHmZcGgzHNxc6nWL0xai5OdIgqBRsvO4ouNo57VxR0jsaF2mTY3wJqm42X29w2ObAu18H44xdVBYlNOLpEXvk7cX5MLBTe6y0Xk9sMkZtGvp2ycfGxkPca2cUQhjW514H9fK+QbaIN4ai9ZmeUfgyA7I3OEkCUd2wOGtJ5KH3AMnH9+gsZM8tB4ICb9wkobGqU4Cao8uTAB506biXlV9r4biMab+CGvgVEc36/bzfQU5zvPwg5vgwEbnD8zh7XBggzMk+bGck48PbeD0QolJgujmEN0UGjZzXqObE3F0t9NALzSyRj6acRUdg5z9zlgOOfuc18ydcCTDWTJ3OAlEQdbJ50mQ8+/ZqC10vsh5jW/rvDZKtporc9rypk3FHTiPPowxNSUs6sQkaqWpOpOqHd7mfqvNcBKOI9ucP1j710P27hPdYHFaWvPTbc4ARdHNoGETiGrs1IA0SPBY4k9+DW1gNSCejh11yr70knsQcvbReeMy2Po0ZO+FnL0/r2EoEZXoJIAJHZzuzDFJThfNmJbOa8Nm1tvC1EoBaVMhIvHu9ZKBLcBVqnqojOPOBf4JBAOvqurj7vYrgYeAzsAAVV3gcc59wC1AEfBrVf3a3d4XeB2IBKYBd9msq6ZWEnH/6Mc7Y2aUpbgYjh50JnPK2s2aBel0Sopzviln7XL+6O1e4fzRO3qIk9p2eAqJcJ7Fh0e7S8yJ14iYk7eHNXSq3kMiS71GODUknq/+TlSKi6AwH4rynem1C/OcRouF+c5yLAfys50aoYIs5zU/25kxsyD7xHp+llOOJclDSe+esoQ2ICY4BsJaQ0J7Z56Zhk2dBKJhUzeRS3RqkuwRhamjAtWm4l7gO1V9XETuddf/6HmAiAQDzwOjgAxgvohMVdVVwArgMuDfpc7pAlwDdAVaANNFpKOqFgEvAuOBuThJxbnAl9X4DMacvoKCnJqIqMbQrBu7d4TQaXha2ccWF8HRw84fz9wDHou7nnfE+eOanwX5mU41fn6W03gwP5NyE5IK4wt15oEICgEJPvH++LagE+ta7NTOaPGJ93iuu0tx4YmkQYuqVm6hDZzkKCzKmTUzPMZppxDZqPIlvCE/paeTllZOORtTD3gzS2lbP9x3DJDmvp8MpFMqqQAGABtUdROAiLzrnrdKVVe728q67ruqmg9sFpENwAAR2QLEqOoc97w3gEuwpMIY5w93VIKzkHJq56o63+zzs5xv9oVH4Vhe5a/Fhc5SkgwUFzrJzUnbitzkQJwk4/gipV6DnGOCQpyZZ0PCnQnjQsJKvUaceB8W5SYO0SeSiLAom/DKmGoqN6kQkT+o6pPu+ytV9QOPff+nqvdX475NVXUXgKruEpEmZRyTBGz3WM8ABpZxXOlz5pY6Jwk45r4vvb1MIjIep1aDxMRE0tPTK7mtqY7s7GwrYz8LTBlHuEspQe7ib4Xukl/WziIg0118x36W/c/K+PRWUU3FNcCT7vv7gA889p0LVJhUiMh0oKzZl/7kZWxlPXStrJ61vHNO6Vqq+jLwMkBqaqpadaZ/pVuVsd9ZGdcMK2f/szI+vVWUVEg578ta/xlVPbvcC4vsEZHmbi1Fc2BvGYdlAK081lsCOyu5bXnnZLjvT+VaxhhjjDkFFVVCajnvy1o/VVOBce77ccBnZRwzH0gRkbYiEoZTc1LZUOFTgWtEJFxE2uI8IJ7nPmrJEpFB4jTEuKGcexpjjDGmiiqqqegpIpk4tRKR7nvc9er2h3oceF9EbgG2AVcCiEgLnK6j56tqoYhMBL7G6VI6SVVXusddCjwHJAJfiMgSVT1HVVeKyPvAKpynqXe4PT8AfsWJLqVfYo00jTHGGJ8qN6lQVb81g1bVA8DIMrbvBM73WJ+G0/2z9HGfAJ+Uc+3HgMfK2L4AKGPoQmOMMcb4gk1PZ4wxxhifsKTCGGOMMT5hSYUxxhhjfEJs+ouKiUgWsDbQcdRxjYH9gQ6ijrMyrhlWzv5nZex/qaoaXZUTbRq8yq1V1X6BDqIuE5EFVsb+ZWVcM6yc/c/K2P9EZEHlR5XNHn8YY4wxxicsqTDGGGOMT1hSUbmXAx1APWBl7H9WxjXDytn/rIz9r8plbA01jTHGGOMTVlNhjDHGGJ+wpMIYY4wxPmFJBSAi54rIWhHZICL3lrFfRORZd/8yEekTiDhrOy/K+Tq3fJeJyGwR6RmIOGuzysrY47j+IlIkIlfUZHx1gTdlLCJpIrJERFaKyA81HWNt58XvilgR+VxElrplfFMg4qzNRGSSiOwVkRXl7K/a3z1VrdcLzgyoG4F2QBiwFOhS6pjzcWY1FWAQ8FOg465ti5flPARo5L4/z8rZ92Xscdz/cCbruyLQcdemxcuf4zicmZJbu+tNAh13bVq8LOP7gSfc94nAQSAs0LHXpgUYDvQBVpSzv0p/96ymAgYAG1R1k6oWAO8CY0odMwZ4Qx1zgTgRaV7TgdZylZazqs5W1UPu6lygZQ3HWNt587MMcCfwEbC3JoOrI7wp42uBj1V1G4CqWjmfGm/KWIFoERGgIU5SUVizYdZuqjoDp9zKU6W/e5ZUQBKw3WM9w912qseYip1qGd6CkyUb71VaxiKSBFwKvFSDcdUl3vwcdwQaiUi6iCwUkRtqLLq6wZsy/hfQGdgJLAfuUtXimgmv3qjS3z0bptup2imtdD9bb44xFfO6DEVkBE5SMdSvEdU93pTxM8AfVbXI+ZJnTpE3ZRwC9AVGApHAHBGZq6rr/B1cHeFNGZ8DLAHOAtoD34rITFXN9HNs9UmV/u5ZUuFkX6081lviZL+neoypmFdlKCI9gFeB81T1QA3FVld4U8b9gHfdhKIxcL6IFKrqpzUSYe3n7e+L/aqaA+SIyAygJ2BJhXe8KeObgMfVefi/QUQ2A52AeTUTYr1Qpb979vgD5gMpItJWRMKAa4CppY6ZCtzgtoYdBBxR1V01HWgtV2k5i0hr4GPgevtWVyWVlrGqtlXVZFVNBj4EbreE4pR48/viM2CYiISISANgILC6huOszbwp4204NUGISFMgFdhUo1HWfVX6u1fvaypUtVBEJgJf47Q6nqSqK0Vkgrv/JZxW8ucDG4BcnCzZnAIvy/kBIAF4wf0mXag2G6HXvCxjUw3elLGqrhaRr4BlQDHwqqqW2W3P/JyXP8d/AV4XkeU41fR/VFWbDv0UiMg7QBrQWEQygAeBUKje3z0bptsYY4wxPmGPP4wxxhjjE5ZUGGOMMcYnLKkwxhhjjE9YUmGMMcYYn7CkwhhjjDE+YUmFMfWcO1vpEo+l3NlNT/G600Qk7lT3VXLNG0XkX9WNzRjjH/V+nApjDEdVtZevL6qq55fe5k4AJWXtM8bUflZTYYwpk4hsEZH/E5E5IrJARPqIyNcisrFkICIRSRORGSLyiYisEpGXRCTI4/zGIpIsIqtF5AVgEdCqZJ973A0iskxElorIm+62i0TkJxFZLCLT3VETK4r1IRGZLCLfuNe+TESeFJHlIvKViIS6x/UVkR/cib6+Lpl1UUR+KSLz3Rg+ckfCREReF5FnRWS2iGwSkSv8Vd7G1AWWVBhjIks9/rjaY992VR0MzAReB64ABgGPeBwzAPgt0B1ncqfLyrhHKs40yr1VdWvJRhHpCvwJOEtVewJ3ubtmAYNUtTfO1Nd/8OJztAcuwJmy+S3ge1XtDhwFLnATi+eAK1S1LzAJeMw992NV7e/GsBpnQrsSzXEmt7sQeNyLOIypt+zxhzGmoscfJXMuLAcaqmoWkCUieR5tIuap6iY4PvTvUJx5RTxtVdW5ZVz/LODDkiGWVfWgu70l8J5bkxAGbPbic3ypqsfcoZuDga88Yk/GSWy64cxoiXtMyVwG3UTkUSAOaIgzRHSJT91ptVdVVmNiTH1nSYUxpiL57muxx/uS9ZLfH6XH+i9r7P+ccq4v5Rz/HPC0qk4VkTTgIW9jVdViETmmJ+YgKIlVgJVuzUtprwOXqOpSEbkRZ06Ek67rEa8xphz2+MMYU10D3Bklg4CrcR5deOs74CoRSQAQkXh3eyyww30/zkdxrgUSRWSwe69Q9/ELQDSwy31Ecp2P7mdMvWNJhTGmdJuKU203MAenrcEKnMcUn3h7oqquxGnX8IOILAWednc9BHwgIjMBn8w+qaoFOG1CnnDvtQQY4u7+M/AT8C2wxhf3M6Y+sllKjTFV5j6a+J2qXhjgUIwxpwGrqTDGGGOMT1hNhTHGGGN8wmoqjDHGGOMTllQYY4wxxicsqTDGGGOMT1hSYYwxxhifsKTCGGOMMT7x/wF1DZgGw7lvngAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "upper_margins = np.array(uppers) - np.array(averages)\n", "lower_margins = np.array(lowers) - np.array(averages)\n", "plt.subplots(figsize=(8, 4))\n", "plt.plot(averages, upper_margins, label='Upper margin of error')\n", "plt.plot(averages, lower_margins, label='Lower margin of error')\n", "plt.xlabel('Empirical mean')\n", "plt.ylabel('Error on the parameter of the distribution')\n", "plt.xlim(0, 1)\n", "plt.grid(True)\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some examples:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "ExecuteTime": { "end_time": "2021-02-16T07:46:53.251662Z", "start_time": "2021-02-16T07:46:53.171875Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Lower boundEmpirical meanUpper boundWidth of the confidence interval
k
00.0000%0.0000%0.0300%0.0300%
10.0003%0.0100%0.0557%0.0555%
20.0024%0.0200%0.0722%0.0698%
30.0062%0.0300%0.0876%0.0815%
40.0109%0.0400%0.1024%0.0915%
50.0162%0.0500%0.1166%0.1004%
60.0220%0.0600%0.1305%0.1085%
70.0281%0.0700%0.1442%0.1160%
80.0345%0.0800%0.1576%0.1230%
90.0412%0.0900%0.1708%0.1296%
100.0480%0.1000%0.1838%0.1359%
200.1222%0.2000%0.3087%0.1865%
300.2025%0.3000%0.4280%0.2255%
400.2859%0.4000%0.5443%0.2584%
500.3713%0.5000%0.6587%0.2873%
600.4582%0.6000%0.7717%0.3135%
700.5461%0.7000%0.8836%0.3375%
800.6348%0.8000%0.9947%0.3598%
900.7243%0.9000%1.1051%0.3808%
1000.8144%1.0000%1.2150%0.4006%
2001.7347%2.0000%2.2938%0.5591%
3002.6744%3.0000%3.3533%0.6789%
4003.6244%4.0000%4.4027%0.7783%
10009.4187%10.0000%10.6047%1.1860%
200019.2198%20.0000%20.7977%1.5778%
300029.1028%30.0000%30.9089%1.8061%
400039.0378%40.0000%40.9680%1.9301%
500049.0151%50.0000%50.9849%1.9697%
600059.0320%60.0000%60.9622%1.9301%
700069.0911%70.0000%70.8972%1.8061%
800079.2023%80.0000%80.7802%1.5778%
900089.3953%90.0000%90.5813%1.1860%
910090.4221%91.0000%91.5539%1.1318%
920091.4510%92.0000%92.5244%1.0735%
930092.4823%93.0000%93.4925%1.0102%
940093.5166%94.0000%94.4576%0.9410%
950094.5545%95.0000%95.4190%0.8645%
960095.5973%96.0000%96.3756%0.7783%
970096.6467%97.0000%97.3256%0.6789%
980097.7062%98.0000%98.2653%0.5591%
990098.7850%99.0000%99.1856%0.4006%
991098.8949%99.1000%99.2757%0.3808%
992099.0053%99.2000%99.3652%0.3598%
993099.1164%99.3000%99.4539%0.3375%
994099.2283%99.4000%99.5418%0.3135%
995099.3413%99.5000%99.6287%0.2873%
996099.4557%99.6000%99.7141%0.2584%
997099.5720%99.7000%99.7975%0.2255%
998099.6913%99.8000%99.8778%0.1865%
999099.8162%99.9000%99.9520%0.1359%
999199.8292%99.9100%99.9588%0.1296%
999299.8424%99.9200%99.9655%0.1230%
999399.8558%99.9300%99.9719%0.1160%
999499.8695%99.9400%99.9780%0.1085%
999599.8834%99.9500%99.9838%0.1004%
999699.8976%99.9600%99.9891%0.0915%
999799.9124%99.9700%99.9938%0.0815%
999899.9278%99.9800%99.9976%0.0698%
999999.9443%99.9900%99.9997%0.0555%
1000099.9700%100.0000%100.0000%0.0300%
\n", "
" ], "text/plain": [ " Lower bound Empirical mean Upper bound Width of the confidence interval\n", "k \n", "0 0.0000% 0.0000% 0.0300% 0.0300%\n", "1 0.0003% 0.0100% 0.0557% 0.0555%\n", "2 0.0024% 0.0200% 0.0722% 0.0698%\n", "3 0.0062% 0.0300% 0.0876% 0.0815%\n", "4 0.0109% 0.0400% 0.1024% 0.0915%\n", "5 0.0162% 0.0500% 0.1166% 0.1004%\n", "6 0.0220% 0.0600% 0.1305% 0.1085%\n", "7 0.0281% 0.0700% 0.1442% 0.1160%\n", "8 0.0345% 0.0800% 0.1576% 0.1230%\n", "9 0.0412% 0.0900% 0.1708% 0.1296%\n", "10 0.0480% 0.1000% 0.1838% 0.1359%\n", "20 0.1222% 0.2000% 0.3087% 0.1865%\n", "30 0.2025% 0.3000% 0.4280% 0.2255%\n", "40 0.2859% 0.4000% 0.5443% 0.2584%\n", "50 0.3713% 0.5000% 0.6587% 0.2873%\n", "60 0.4582% 0.6000% 0.7717% 0.3135%\n", "70 0.5461% 0.7000% 0.8836% 0.3375%\n", "80 0.6348% 0.8000% 0.9947% 0.3598%\n", "90 0.7243% 0.9000% 1.1051% 0.3808%\n", "100 0.8144% 1.0000% 1.2150% 0.4006%\n", "200 1.7347% 2.0000% 2.2938% 0.5591%\n", "300 2.6744% 3.0000% 3.3533% 0.6789%\n", "400 3.6244% 4.0000% 4.4027% 0.7783%\n", "1000 9.4187% 10.0000% 10.6047% 1.1860%\n", "2000 19.2198% 20.0000% 20.7977% 1.5778%\n", "3000 29.1028% 30.0000% 30.9089% 1.8061%\n", "4000 39.0378% 40.0000% 40.9680% 1.9301%\n", "5000 49.0151% 50.0000% 50.9849% 1.9697%\n", "6000 59.0320% 60.0000% 60.9622% 1.9301%\n", "7000 69.0911% 70.0000% 70.8972% 1.8061%\n", "8000 79.2023% 80.0000% 80.7802% 1.5778%\n", "9000 89.3953% 90.0000% 90.5813% 1.1860%\n", "9100 90.4221% 91.0000% 91.5539% 1.1318%\n", "9200 91.4510% 92.0000% 92.5244% 1.0735%\n", "9300 92.4823% 93.0000% 93.4925% 1.0102%\n", "9400 93.5166% 94.0000% 94.4576% 0.9410%\n", "9500 94.5545% 95.0000% 95.4190% 0.8645%\n", "9600 95.5973% 96.0000% 96.3756% 0.7783%\n", "9700 96.6467% 97.0000% 97.3256% 0.6789%\n", "9800 97.7062% 98.0000% 98.2653% 0.5591%\n", "9900 98.7850% 99.0000% 99.1856% 0.4006%\n", "9910 98.8949% 99.1000% 99.2757% 0.3808%\n", "9920 99.0053% 99.2000% 99.3652% 0.3598%\n", "9930 99.1164% 99.3000% 99.4539% 0.3375%\n", "9940 99.2283% 99.4000% 99.5418% 0.3135%\n", "9950 99.3413% 99.5000% 99.6287% 0.2873%\n", "9960 99.4557% 99.6000% 99.7141% 0.2584%\n", "9970 99.5720% 99.7000% 99.7975% 0.2255%\n", "9980 99.6913% 99.8000% 99.8778% 0.1865%\n", "9990 99.8162% 99.9000% 99.9520% 0.1359%\n", "9991 99.8292% 99.9100% 99.9588% 0.1296%\n", "9992 99.8424% 99.9200% 99.9655% 0.1230%\n", "9993 99.8558% 99.9300% 99.9719% 0.1160%\n", "9994 99.8695% 99.9400% 99.9780% 0.1085%\n", "9995 99.8834% 99.9500% 99.9838% 0.1004%\n", "9996 99.8976% 99.9600% 99.9891% 0.0915%\n", "9997 99.9124% 99.9700% 99.9938% 0.0815%\n", "9998 99.9278% 99.9800% 99.9976% 0.0698%\n", "9999 99.9443% 99.9900% 99.9997% 0.0555%\n", "10000 99.9700% 100.0000% 100.0000% 0.0300%" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ks_examples = (\n", " list(range(10)) \n", " + list(range(10, 100, 10)) \n", " + list(range(100, 500, 100))\n", " + list(range(1000, 9000, 1000))\n", " + list(range(9000, 9900, 100))\n", " + list(range(9900, 9990, 10))\n", " + list(range(9990, 10001, 1))\n", ")\n", "df = pd.DataFrame()\n", "df.index.name = 'k'\n", "for k in ks_examples:\n", " df.loc[k, 'Lower bound'] = \"{:.4%}\".format(lowers[k])\n", " df.loc[k, 'Empirical mean'] = \"{:.4%}\".format(averages[k])\n", " df.loc[k, 'Upper bound'] = \"{:.4%}\".format(uppers[k])\n", " df.loc[k, 'Width of the confidence interval'] = \"{:.4%}\".format(uppers[k] - lowers[k])\n", "df" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }